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We perform a quantitative simulation of the repulsive Fermi-Hubbard model using an ultra-cold 
gas trapped in an optical lattice. The entropy of the system is determined by comparing accurate 
measurements of the equilibrium double occupancy with theoretical calculations over a wide range 
of parameters. We demonstrate the applicability of both high-temperature series and dynamical 
mean-field theory to obtain quantitative agreement with the experimental data. The reliability of 
the entropy determination is confirmed by a comprehensive analysis of all systematic errors. In the 
center of the Mott insulating cloud we obtain an entropy per atom as low as 0.77/cb which is about 
twice as large as the entropy at the Neel transition. The corresponding temperature depends on the 
atom number and for small fillings reaches values on the order of the tunneling energy. 

PACS numbers: 05.30.Fk, 03.75.Ss, 67.85.-d, 71.10.Fd 



Experimental progress in the field of atomic quantum 
gases has led to a new approach to quantum many-body 
physics. In particular, the combination of quantum de- 
generate and strongly interacting Fermi gases [TJ [5] with 
optically induced lattice potentials [Hj now allows the 
study of a centerpiece of quantum condensed matter 
physics, the Fermi-Hubbard model [4]. The high level 
of control over the atomic systems has led to the con- 
cept of quantum simulation, which for the case of the 
Fermi-Hubbard model is expected to provide answers to 
intriguing open questions of frustrated magnetism and 
d-wave superfluidity [5] . Recent experiments [6j [7] have 
indeed demonstrated that the strongly correlated regime 
of the repulsive Fermi-Hubbard model is experimentally 
accessible and the emergence of a Mott insulating state 
has been observed. In this Letter, we succeed in per- 
forming a quantitative simulation of the Fermi-Hubbard 
model using cold atoms. The level of precision of the 
experiment enables us to determine the entropy and the 
temperature of the system, and thereby to quantify the 
approach to the low temperature phases. 

The main challenge for the quantum simulation of 
the Fermi-Hubbard model is a further reduction in tem- 
perature. Here the lack of a quantitative thermometry 
method in the lattice is a key obstacle. For strongly 
correlated bosonic systems thermometry has recently 
been demonstrated by direct comparison with quantum 
Monte-Carlo simulations [5] or by using the boundary 
of two spin polarized clouds [Sj. In the fermionic case, 
previous methods to determine the temperature could 
be used only in limiting regimes of the Hubbard model, 
namely the noninteracting |10l 111] and zero-tunneling 



[BJ [12] regimes. However, intermediate interactions are 
most interesting for quantum simulation of the Fermi- 
Hubbard model and no reliable thermometry method has 
been available up to now. 

In both the metallic and Mott insulating regimes an 
accurate measurement of the double occupancy provides 
direct access to thermal excitations. We analyze the 
crossover from thermal creation of double occupancies to 
thermal depletion which is unique to a trapped system 
(see Fig. [T]) . The variability of the double occupancy with 
respect to temperature allows the entropy of the system 
to be inferred directly by comparison with two ab-initio 
theoretical methods. By determining all other quantities 
entering the analysis separately and with methods that 
are independent of the double occupancy measurement, 
we demonstrate the versatility of the double occupancy 
in quantifying the state of the system. 

To obtain a quantum degenerate Fermi gas we adhere 
to the procedure described in previous work [BJ. A bal- 
anced spin mixture of K atoms in the mp — —9/2 and 
—5/2 magnetic sublevels of the F = 9/2 hyperfine man- 
ifold is evaporatively cooled in a crossed beam optical 
dipole trap, with less than 1.2% of the atoms remaining 
in the mp — — 7/2 state. We prepare Fermi gases with 
total atom numbers between N = 30 x 10 3 and 300 x 10 3 . 
The atom number is calibrated using strong saturation 
imaging |13j at high magnetic field, with a systematic 
error ~ 10%. 

The optical lattice potential is then ramped up in 0.2 s 
and has a simple cubic symmetry with lattice constant 
d = 532 nm. Its depth is determined from Raman-Nath 
diffraction of 87 Rb and confirmed by resonant excitation 
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FIG. 1: Double occupancy: experiment versus theory. Points 
and error bars are the mean and standard deviation of at 
least three experimental runs. The solid curve in each panel 
is the best fit of the second order high-temperature series 
to the experimental data and yields specific entropies of 
s = 2.2(2), 2.0(5), 1.9(4), 1.6(4) for the different interactions 
strengths of U/6t = 1.4(2), 2.4(4), 3.2(5), 4.1(7). Curves for 
s — 1.3 (dashed curve) and 2.5 (dotted curve) represent the in- 
terval of specific entropy measured before and after the ramp- 
ing of the lattice. We use fee = 1- 



of atoms to higher bands [H]. In the lowest Bloch band 
the tunneling matrix element is t/h = 174(30) Hz where 
h is Planck's constant. From transverse oscillations in 
the standing wave potentials of each lattice beam and the 
dipole trap we extract overall harmonic trapping frequen- 
cies of u) x , y , z /2ir = [49.4(9), 52.6(6), 133.0(10)] Hz and a 
geometric mean of uj/2tt = 70.1(5) Hz. The characteris- 
tic filling is p = N/N [15] where N = (12t/muj 2 d 2 ) 3 / 2 
and m is the 40 K atom mass. The characteristic atom 
number No relates the bandwidth to the inhomogeneity 
and corresponds to the atom number per spin that yields 
half filling in the center of the trap at zero temperature 
without interaction. 

We tune the interaction between the atoms by adjust- 
ing the scattering length in the proximity of an s-wave 
Feshbach resonance before loading into the lattice. To 
determine the scattering length, the width of the reso- 
nance was measured by using the suppressed dephasing 
of Bloch oscillations [TS] to locate the zero crossing of the 
scattering length. We obtain a width of the resonance of 
7.5(1) G which deviates from the result of Ref. [17] where 
the mean-field energy was measured. We infer the on-site 
interaction energy U from the scattering length and the 
Wannier function in the lowest Bloch band [IS]. This ab 
initio U is experimentally validated using resonant exci- 
tation of double occupancies by lattice modulation [6lll9|. 
We cover the range from weak repulsion in the metallic 
regime to strong repulsion with a Mott insulating core 
using scattering lengths between 200eto and 650ao, where 
ao is the Bohr radius. We choose values of the Hubbard 



parameter U/6t = 1.4(2), 2.4(4), 3.2(5) and 4.1(7). Be- 
cause of the lattice loading process, beam intensity noise 
and incoherent photon scattering, the atoms heat up dur- 
ing preparation. Before loading into the lattice, the tem- 
perature in the dipole trap is around 0.13Tp independent 
of the atom number as determined from the momentum 
distribution after time-of-flight. Here Tp is the Fermi 
temperature. This corresponds to an entropy per atom 
of s = S/N w 1.3 [201 • Since the system is isolated from 
the environment, the temperature changes significantly 
even when adiabatically loading into the lattice. The 
entropy, however only changes due to nonadiabatic pro- 
cesses. Therefore we can find a typical upper limit of the 
specific entropy in the lattice by reversing the loading 
sequence and subsequently measuring the temperature. 
Here we obtain s < 2.5. 

After loading the atoms into the lattice we determine 
the double occupancy. A sudden increase of the lat- 
tice depth suppresses further tunneling. The fraction of 
atoms on doubly occupied lattice sites D is then obtained 
by combining rf spectroscopy, Stern-Gerlach separation 
of the spin components and absorption imaging [H] I19j . 
Here we account for the independently determined offset 
due to the imperfection of the initial spin mixture. From 
long term reproducibility and comparison with the adi- 
abatic formation of molecules via magnetic field sweeps 
we conclude that the relative systematic uncertainty of 
the double occupancy measurement is 10%. 

Because of the harmonic trapping potential, the tem- 
perature behavior of the double occupancy can be 
markedly different from that of homogeneous systems 
[10[ l2Tj . In a homogeneous system the double occupancy 
increases with temperature in most regimes of filling and 
interaction strength due to thermal activation. However, 
in a harmonically trapped system an increase in temper- 
ature allows the atoms to reach outer regions of the trap, 
in turn reducing the density in the central region: in this 
case thermal excitations do not populate doubly occupied 
states but rather deplete them through the decrease of 
the density. The regimes depicted in Fig. [T] demonstrate 
the competition between thermal activation and the ef- 
fect of the trapping potential on the double occupancy as 
a function of filling and entropy. To extract the entropy 
the experimental data are compared with theoretical re- 
sults. The curves in Fig. [T] correspond to the best fitting 
entropy and its experimental bounds. 

We apply a high-temperature series expansion [22] as 
well as single-site dynamical mean-field theory (DMFT) 
[23] with a continuous time quantum Monte Carlo solver 
[24) . In the experimentally relevant regime we find the 
high-temperature series and DMFT to be in agreement 
to within 0.2%. For simplicity, the theoretical curves 
shown in this Letter are therefore generated using the 
second order high-temperature series unless stated oth- 
erwise. The entropy is determined from a one-parameter 
least-squares fit of the high-temperature series D(s,pi) 
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to the experimental data points Di weighting them ac- 
cording to their statistical errors <tdi- The fit minimizes 
X 2 = Yli{D{s, pi) — Di) 2 ja 2 Di . The series is able to ac- 
curately reproduce the measured double occupancy for 
all shown interaction strengths. We find deviations of 
the experimental data only at the lowest fillings at low 
repulsion, indicating that for very small atom numbers 
and weak interaction the cooling and loading procedure 
may fail to produce a constant entropy per atom. 

The size and direction of the corridors between the 
initial and the final entropy in the dipole trap has impli- 
cations for the usefulness of the double occupancy when 
performing thermometry. In Fig. [lja) the behavior for 
low repulsion U/6t = 1.4 is shown. With increasing fill- 
ing the system transforms from a dilute gas to an in- 
creasingly dense metal with high double occupancy. In 
this case the effect of the trapping potential dominates 
and D decreases with increasing entropy. Because of its 
large \8D/ds\, the regime of Fig. [lja) is well suited for 
thermometry 

At intermediate repulsion strengths in Fig. [ljb) and 
(c), double occupancies become increasingly suppressed 
and \dD/ds\ decreases. For each interaction strength 
dD/ds changes sign at a certain filling. These points 
mark the crossover to thermal suppression of double 
occupancies. If dD/ds vanishes, the theory becomes 
parameter-free to first order and can be used to further 
determine other calibration factors, e.g., the characteris- 
tic filling. 

Fig. [ijd) shows data for clouds in the Mott insulating 
regime. It exhibits a pronounced plateau of suppressed 
double occupancy at intermediate fillings owing to a van- 
ishing core compressibility, a characteristic signature of 
a Mott insulating core [5J [T3] . Large filling can increase 
the chemical potential to values comparable with U and 
thus create double occupancy. In this regime the ther- 
mal activation of double occupancies dominates over the 
thermal decrease of density due to the trapping potential. 
If dD/ds > a large fraction of particles resides in the 
Mott insulating core. Here the chemical potential is high 
enough to prevent holes from entering the center and ad- 
ditionally the density of states is sufficiently gapped to 
allow only few thermally excited double occupancies. 

We now consider the systematic errors of all param- 
eters and measurements to assess the absolute reliabil- 
ity of the present method in determining the entropy. 
Table [I] lists the contributions. The sensitivity of the 
least-squares fit to variation of the respective parameter 
shows the sign of the influence as well as the magnitude. 
The total relative uncertainties are below 25% for all four 
interaction strengths which confirms the validity of the 
determined entropies. It is apparent that the systematic 
errors dominate and that especially the atom number and 
double occupancy calibrations are critical. The observed 
increase of the specific entropy with decreasing interac- 
tion can be explained by an interaction-dependent global 
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TABLE I: Error budget of the entropy determination. The 
table lists the sensitivity of the fit dsRt/d(-) to the changes 
in the system's parameters scaled by their systematic errors 
Sr.). For a positive contribution an increase in the parameter 
would lead to an increase in the apparent entropy. The con- 
tributions are added in quadrature to the fit error estimate 
a'i = 2(d 2 x 2 /ds 2 )' 1 to obtain the total uncertainty of the 
entropy. 



adiabaticity of the preparation [25] or by a combination 
of systematic errors in N and D. 
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FIG. 2: Properties of experimental regimes and validity of 
theoretical methods. Panels (a)-(c) show the central density 
no, central compressibility kq and temperature T reached in 
the corresponding Hubbard model as a function of characteris- 
tic filling p for the parameters of Fig.JT] Panel (d): Agreement 
between high-temperature series (HTS, second order dotted 
line, sixth order dash-dotted, tenth order dashed line) and 
DMFT (solid line) for U /U = 2.5 and p = 16.5 as a func- 
tion of temperature in the lattice T/6t. For low temperatures 
T <t the series starts to diverge. 

From the theoretical description, several unique prop- 
erties of trapped repulsively interacting Fermi-Hubbard 
systems can be derived. Figures |2|a) and (b) show the 
central density no and compressibility kq — dng/dp ver- 
sus characteristic filling for the interaction strengths and 
specific entropies of Fig. [I] The plateau in no and the 
reduction of ko f° r V /6t = 4.1 are signatures of the Mott 
insulating regime [5T] . Compared to the result for a non- 
interacting system where the ground state has a com- 
pressibility of 1.69/6i at half filling, the compressibil- 
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ity is suppressed by a factor of 50 to values as low as 
k « 0.03/6i 

The entropy as determined above needs to be related 
to a temperature to allow for comparison with models 
of homogeneous systems. Figure [2jc) shows this tem- 
perature in units of the half bandwidth as a function of 
characteristic filling. The behavior is similar to that of a 
Fermi gas in a harmonic trap where the temperature at 
constant specific entropy increases with the atom number 
[TUl ED] . At the lowest fillings of p — 5 the temperature 
in the lattice even approaches the energy scale of the 
tunneling T ~ t. At these low temperatures the results 
of high-temperature series and DMFT start to deviate 
considerably, see Fig. [2jd) . 




s i H ln4h s i 



reaches values close to the maximum of Si = In 4 in some 
regions. For large repulsion the sites in the perimeter 
of the cloud where rii ~ 2/3 carry most of the entropy. 
This can be understood in the atomic limit at large U . 
At Hi — 2/3 each site has three equally likely states and 
can accommodate Sj — In 3 of entropy. The Mott insulat- 
ing core can only absorb In 2 of spin entropy. Mean-field 
theory of the Heisenberg model predicts this to coincide 
with the entropy at its critical point. However, quan- 
tum fluctuations lower the entropy at the Neel temper- 
ature where magnetic long-range order sets in to about 
sh ~ In 2/2 sa 0.35 [2B]. We have verified numerically 
that for the Heisenberg model with exchange coupling J 
the entropy is s = 0.338(5) at T N6cl /J = 0.946(1) [27l [28] 
by integrating the energy, S — J dE/T. Integration from 
above or below the Neel temperature agree. Additionally, 
we performed a new study of the Hubbard model using a 
diagrammatic determinant Monte Carlo method [29] . For 
U/t = 8, the critical temperature is Tweei/t = 0.325(7) 
[52] , and the critical entropy s N6ol = 0.345(45). This dif- 
fers from the mean-field calculation including fluctuation 
corrections of Ref. [30] and is a factor of 2 less than the 
experimental results presented here. 

We acknowledge funding from SNF, NCCR MaNEP, 
NAME-QUAM (EU, FET-Open), SCALA (EU), ANR 
(FABIOLA and FAMOUS), DARPA-OLE, and "Triangle 
de la Physique" . 



FIG. 3: Density and entropy distribution in the trap. For the 
interaction strengths and entropies of Fig. [I] the density rn 
(upward) and entropy s 4 (downward) per site i at a charac- 
teristic filling of p — 15 in a spherically symmetric system are 
shown. The buffering effect of the low-density shell around the 
Mott insulating core becomes clearly visible for U/Qt = 4.1. 
There, the entropy reaches values of about twice the critical 
entropy of the Heisenberg model sh ~ In 2/2. 

In a finite trapped system, antiferromagnetic order re- 
quires two conditions to be met. It can only be estab- 
lished in regions of sufficiently constant atom density and 
low specific entropy. The experimental situation with 
respect to those conditions is depicted in Fig. [3] for a 
characteristic filling of p = 15 and the same parameters 
as those shown in Fig. [T] The upward axes show the 
spatial density distribution. At low repulsive interaction 
[Figs.[3]ja) and[3]jb)] the system has a density above one 
in the center which corresponds to a significant doping. 
In Fig.^c), double occupancies and holes compensate in 
the center and lead to an average density close to unity 
which then starts to deviate a few sites away from the 
center. Only for the largest repulsion Fig.|3jd), the Mott 
insulating core is robust against the confining potential 
over the central 20 sites. Here, the density changes by 
less than 1%. 

The entropy per site is shown on the downward axes in 
Fig. [3] It is highest in regions of the trap where the den- 
sity is most variable. For small interaction strengths it 
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